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Abstract 

This paper summarizes the results obtained with the 
NSU3D unstructured multigrid solver for the AIAA 
Drag Prediction Workshop held in Anaheim. CA, June 
2001. The test case for the workshop consists of a 
wing-body configuration at transonic flow conditions. 
Flow analyses for a complete test matrix of lift coeffi- 
cient values and Mach numbers at a constant Reynolds 
number are performed, thus producing a set of drag 
polars and drag rise curves which are compared with 
experimental data. Results were obtained indepen- 
dently by both authors using an identical baseline grid, 
and different refined grids. Most cases were run in 
parallel on commodity cluster- type machines while the 
largest cases were run on an SGI Origin machine using 
128 processors. The objective of this paper is to study 
the accuracy of the subject unstructured grid solver 
for predicting drag in the transonic cruise regime, to 
assess the efficiency of the method in terms of conver- 
gence. epu time' and memory, and to determine the 
effects of grid resolution on this predictive' ability and 
its computational efficiency. A good predictive abil- 
ity is demonstrated over a wide range of conditions, 
although accuracy was found to degrade for cases at 
higher Mach numbers and lift value's where increas- 
ing amounts of flow separation occur. The ability to 
rapidly compute large numbers of cases at varying flow 
conditions using an unstructured solver on inexpensive 
clusters of commodity computers is also demonstrated. 

Introduction 

Computational fluid dynamics has progressed to the 
point where Reynolds- averaged Navier-Stokes solvers 
have become standard simulation tools for predict- 
ing aircraft aerodynamics. These solvers are routinely 
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used to predict aircraft force coefficients such as lift, 
drag and moments, as well as the changes in these 
values with design changes. In order to be useful to 
an aircraft designer, it is generally acknowledged that 
the computational method should be capable of pre- 
dicting drag to within several counts. While Reynolds- 
averaged Navier-Stokes solvers have made great strides 
in accuracy and affordability over the last decade, the 
stringent accuracy requirements of the drag prediction 
task have proved difficult to achieve. This difficulty is 
compounded by the multitude' of Navier-Stokes solver 
formulations available, as well as by the effects on 
accuracy of turbulence modeling and grid resolution. 
Therefore, a particular Navier-Stokes solver must un- 
dergo extensive validation including the determination 
of adequate grid resolution distribution, prior to being 
trusted as a useful predictive tool. With these issues in 
mind, the AIAA Applied Aerodynamics technical com- 
mittee organized a Drag Prediction Workshop, held 
in Anaheim CA. June 2001. 1 in order to assess the 
predictive capabilities of a number of state-of-the-art 
computational fluid dynamics methods. The chosen 
configuration, denoted as DLR-F4 2 and depicted in 
Figure 1. consists of a wing-bodv geometry, which is 
representative* of a modern supercritical swept wing 
transport aircraft. Participants included Reynolds- 
averaged Navier-Stokes formulations based on block- 
structured grids, overset grids, and unstructured grids, 
thus affording an opportunity to compare* these' meth- 
ods on ail f'qual basis in terms of accuracy and effi- 
ciency. A standard mesh was supplied for each type.' 
of methodology, with participants encouraged to pro- 
duce results on additionally refined meshes, in order to 
assess the effects of grid resolution. A Mach number 
versus lift coefficient (C L ) matrix of test cases was de- 
fined. which included mandatory and optional cases. 
The calculations were initially run by the participants 
without knowledge of the experimental data, and a 
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compilation of all workshop results including a statis- 
tical analysis of these results was performed by the 
committee.' 5 



Fig. 1 Definition of Geometry for Wing-Body 

Test Case (taken from Ref. 2 ) 

This paper describes the results obtained for this 
workshop with the unstructured mesh Navior- Stokes 
solver NS13D. 1-6 This solver has been well validated 
and is currently in use in both a research setting and 
an industrial production environment. Results were 
obtained independently by both authors on the' base- 
line workshop grid, and on two refined grids generated 
independently by both authors. All required and op- 
tional cases were run using the baseline grid and one 
refined grid, while the most highly refined grid was 
only run on the mandatory cast’s. The runs were per- 
formed on three different types of parallel machines at 
two different locations. 

Flow Solver Description 

The XSl 3D code solves the Reynolds averaged 
Navier-Stokes equations on unstructured meshes of 
mixed element types which may include tetrahedra, 
pyramids, prisms, and hexahedra. All elements of the 
grid are handled by a single unifying edge-based data- 
structure in the flow solver. 1 

Tetrahedral elements are employed in regions where 
the grid is nearly isotropic, which generally correspond 
to regions of in viscid flow, and prismatic cells are em- 
ployed in regions dost' to the wall, such as in boundary 
layer regions where the grid is highly stretched. Tran- 
sition between prismatic and tetrahedral cell regions 
occurs naturally when only triangular prismatic faces 
are exposed to the tetrahedral region, but requires a 
small number of pyramidal cells (cells formed by 5 ver- 
tices) in cases where quadrilateral prismatic faces are 
exposed. 

Flow variables are stored at the vertices of the mesh, 
and the governing equations are discretized using a 
central difference finite- volume technique with added 
artificial dissipation. The matrix formulation of the 
artificial dissipation is employed, which corresponds 


to an upwind scheme based on a Roe-Riemann solver. 
The thin-layer form of the Navier-Stokes equations is 
employed in ail cases, and the viscous terms are dis- 
cretized to second-order accuracy by finite-difference 
approximation. 1 For mulrigrid calculations, a first- 
order discretization is employed for the convective 
terms on the coarse grid levels. 

The basic time-stepping scheme is a three-stage ex- 
plicit multistage scheme with stage coefficients opti- 
mized for high frequency damping properties/ and 
a CFL number of 1.8. Convergence is accelerated 
by a local block- Jacobi preconditioner in regions of 
isotropic grid cells, which involves inverting a 5 x 5 
matrix for each vertex at each stage.*" 10 In bound- 
ary layer regions, where the grid is highly stretched, 
a line smoother is employed, which involves inverting 
a block tridiagonal along lines constructed in the un- 
structured mesh by grouping together edges normal to 
the grid stretching direction. The line smoothing tech- 
nique has been shown to relieve the numerical stiffness 
associated with high grid anisotropy. 11 

An agglomeration multigrid algorithm 1,1 * is used 
to further enhance convergence to steady-state. In 
this approach, coarse levels are constructed by fus- 
ing together neighboring fine grid control volumes to 
form a smaller number of larger and more complex 
control volumes on the coarse grid. 'Hi is process is 
performed automatically in a pre-processing stage by 
a graph-based algorithm. A multigrid cycle consists 
of performing a time-step on the fine grid of the se- 
quence, transferring the flow solution and residuals to 
the coarser level, performing a time-step on the coarser 
level, and then interpolating the corrections back from 
the coarse level to update the fine grid solution. The 
process is applied recursively to the coarser grids of 
the sequence. 

The single equation turbulence model of Spalart and 
Allmaras 13 is utilized to account for turbulence ef- 
fects. This equation is discretized and solved in a 
manner completely analogous to the flow equations, 
with the exception that the convective terms are only 
discretized to first-order accuracy. 

The unstructured mulrigrid solver is parallelized 
by partitioning the domain using a standard graph 
partitioned 1 and communicating between the various 
grid partitions running on individual processors us- 
ing either the MPI message- passing library. 1 ’ > or the 
OpenMP compiler directives. 10 Since OpenMP gener- 
ally has been advocated for shared memory architec- 
tures. and MPI for distributed memory architectures, 
this dual strategy not only enables the solver to run 
efficiently on both types of memory architectures, but 
can also be used in a hybrid two-level mode, suitable 
for networked clusters of individual shared memory 
multi-processors. 6 For the results presented in the pa- 
per, the solver was run on distributed memory PC 
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dusters and an SGI Origin 2000. using the MPI pro- 
gramming model exclusively. 

Grid Generation 

The baseline grid supplied for the workshop was gen- 
erated using the VGRIDns package. 1 7 This approach 
produces fully tetrahedral meshes, although it is ca- 
pable of generating highly stretched semi-structured 
tetrahedral elements near the wall in the boundary- 
layer region, and employs moderate spanwise stretch- 
ing in order to reduce the total number of points. A 
semi-span geometry was modeled, with the far-held 
boundary located 50 chords away from the origin, re- 
sulting in a total of 1.65 million grid points. 9.7 million 
tetraliedra. and 56.000 wing-body surface* points. The 
chordwise grid spacing at the leading edge was pre- 
scribed as 0.250 mm and 0.500 mm at. the trailing 
edge, using a dimensional mean chord of 142.1 mm. 
The trailing edge is blunt, with a base* thickness of 
0.5 % chord, and the baseline 1 mesh contained 5 grid 
points across the trailing edge. The normal spacing at 
the wall is 0.001 mm, which is designed to produce a 
grid spacing corresponding to y^ = 1 for a Reynolds 
number of 5 million. A stretching rate of 1.2 was pre- 
scribed for the growth of cells in the. 1 normal direction 
near the wall, in order to obtain a minimum of 20 
points in the boundary layer. 

Because the NSU3D solver is optimized to run on 
mixed element meshes, the fully tetrahedral baseline 
mesh is subsequently converted to a mixed element 
mesh by merging the semi-structured tetrahedral lay- 
ers in the boundary layer region into prismatic ele- 
ments. This is done in a pre-processing phase where 
triplets of tetrahedral layers are identified and merged 
into a single prismatic element, using information iden- 
tifying these elements as belonging to the stretched 
viscous layer region as opposed to the isotropic inviscid 
tetrahedral region. The merging operation results in a 
total of 2 million created prismatic elements, while the 
number of tetrahedral cells is reduced to 3.6 million, 
and a total of 10090 pyramidal elements are created to 
merge prismatic elements to tetrahedral elements in 


regions where quadrilateral faces from prismatic ele- 
ments are adjacent to tetrahedral elements. 

A higher resolution mesh was generated by the sec- 
ond author using VGRIDns with smaller spacings in 
the vicinity of the wing root. rip. and trailing edge, 
resulting in a total of 3 million grid points, and 73.000 
wing-bodv surface points. One of the features of this 
refined grid is the use of a total of 1 7 points across the 
wing trailing edge versus 5 for the baseline grid. After 
th(' merging operation, this grid contained a total of 
3.7 million prisms and 6.6 million tetraliedra. 

An additional fine mesh was obtained by the first 
author through global refinement of the baseline work- 
shop mesh. This strategy operates directly on the 
mixed prismatic-tetrahedral mesh, and consists of sub- 
dividing each element into 8 smaller self-similar ele- 
ments. thus producing an 8:1 refinement of the original 
mesh. 18 The final mesh obtained in this manner con- 
tained a total of 13.1 million points with 16 million 
prismatic elements and 28.8 million tetrahedral ele- 
ments, and 9 points across the blunt trailing edge of 
the wing. This approach can rapidly generate very 
large meshes which would otherwise be very time con- 
suming to construct using the original mesh generation 
software. One drawback of the current approach is 
that newly generated surface points do not lie exactly 
on the original surface description of the model geom- 
etry. but rather along a linear interpolation between 
previously existing surface coarse grid points. For a 
single level of refinement, this drawback is not ex- 
pected to have a noticeable effect on the results. An 
interface for re-projecting new surface points onto the 
original surface geometry is currently under consider- 
ation. 

The baseline grid was found to be sufficient to re- 
solve all major flow features. The computed surface 
pressure coefficient on the baseline grid for a Mach 
number of 0.75, Reynolds number of 3 million, and 
Ci - 0.6 is shown in Figure 2. illustrating good reso- 
lution of the upper surface shock. A small region of 
separation is also resolved in the wing root area, as 
shown by the surface streamlines for the same flow 
conditions, in Figure 3. 


Table 1 Grids and Corresponding Run Times 


Grid 

No. Points 

No. Tets 

No. Prisms 

Memory 

Run Time 

Hardware 

Grid 1 

1.63 x 10° 

2 x 10“ 

3.6 x 10“ 

2.8 Gbytes 

2.6 hours 

16 Pentium IV 1.7GHz 

Grid 1 

1.65 x 10 s 

2 x 10° 

3.6 x 10“ 

2.4 Gbytes 

8 hours 

4 DEC Alpha 21264 (G67MHZ) 

Grid 1 

1.65 x 10“ 

2 x 10“ 

3.6 x 10“ 

3.0 Gbytes 

45 min. 

64 SGI Origin 2000 (400MHz) 

Grid 2 

3.0 x 10° 

3.7 x 10“ 

6.6 x 10“ 

4.2 Gbytes 

8 hours 

8 DEC Alpha 21264 (667MHZ) 

Grid 3 

13 x 10° 

16 x 10“ 

28.8 x 10“ 

27 Gbytes 

4 hours 

128 SGI 02000 (400MHz) 
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Figure* 1 depicts the computed i) r values at the break 
section for the same How conditions, indicating values 
well below unity over the entire lower surface and a ma- 
jority of the upper surface. The convergence history 
for this cast 1 is shown in Figure 3. The How is initialized 
as a uniform How at freest ream conditions, and ten sin- 
gle grid cycles (no multigrid) are employed to smooth 
the initialization prior to the initiation of the multi- 
grid iteration procedure. A total residual reduction of 
approximately 3 orders of magnitude' is achieved over 
300 multigrid cycles. Convergence in the lift coeffi- 
cient is obtained in as little as 200 multigrid cycles 
for this case, although all case's are run a minimum of 
300 multigrid cycles as a conservative convergence cri- 
terion. This convergence behavior is representative of 
the majority of cases run. with some of the high Mach 
number and high Cj, cases involving larger regions of 
separation requiring up to 800 to 1000 multigrid cycles. 
A How solution on the baseline grid requires 2.8 Gbytes 
of memory and a total of 2.0 hours of wall clock time 
(for 300 multigrid cycles) on a cluster of commodity 
components using 1G Pentium IV 1.7 GHz processors 
communicating through 100 Mbit Ethernet. This case 
was also run on 4 DEC Alpha processors, requiring 2.4 
Gbytes of memory and 8 hours of wall clock time. This 
case was also benchmarked on G4 processors (400MHz) 
of an SGI Origin 2000. requiring 3 Gbytes of memory 
and 43 minutes of wall clock time. The memory re- 
quirements are independent of the specific hardware 
and arc only a function of the number of partitions 
used in the calculations. 

The cases using the 3 million point grid wen' run 
on a cluster of 8 DEC Alpha processors communicat- 
ing through 100 Mbit Ethernet and required approx- 
imately 8 hours of wall clock time and 1.2 Gbytes of 
memory. 

The 13 million point grid cast's wore run on an SGI 
Origin 2000. using 128 processors and required 4 hours 
of wall clock time and 27 Gbytes of memory. A descrip- 
tion of the three grids employed and the associated 
computational requirements oil various hardware 1 plat- 
forms is give'll in Table' 1. 



Fig. 2 Baseline Grid and Computed Pressure 
Contours at Mach=0.75, Cr. — 0.6. Re = 3 mil- 
lion 
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Fig. 3 Computed Surface Oil Flow Pat- 
tern in Wing Root Area on Baseline Grid for 
Maeh=0.75, Cl = 0.6, Re = 3 million 


Fig. 5 Density Residual and Lift Coefficient 
Convergence History as a Function of Multigrid 
Cycles on Baseline Grid for Mach=0.75, Cl = 
0.6, Re — 3 million 


Lower Surface 
Upper Surface 



0,3 E“ 


i L_J 

0 0.25 0.5 0 75 1 

x/c 


Fig. 4 Computed on wing surface at span 
break section on baseline grid for Mach=0.75, 
Cl — 0.6, Re = 3 million 


Results 

The workshop test eases comprised two required 
cases and two optional cases. These cases are described 
in Table 2. For all cases the Reynolds number is 3 mil- 
lion. The first tost case is a single point at Mat'll 
0.75 and Cj, — 0.5. The second test cast 1 involves the 
computation of the drag polar at Mach =0.75 using in- 
cidences from -3.0 to +2.0 degrees in increments of 1 
degree. Tht 1 optional Cases 3 and 4 involve 1 a matrix 
of Mach and Cl values in order to compute drag rise* 
curves. Since an automated approach for computing 
fixed Cl cases has not been implemented, a complete 
drag polar for each Mach number was computed for 
Cases 3 and 4. For the 1 baseline 1 grid, the incidence for 
the 1 prescribed lift value 1 was them interpellated from 
the drag polar using a cubic spline fit. and the How 
was recomputed at this prescribed incidence. The fi- 
nal force coefficients were then interpolated from the 
value's computed in this ease 1 to the prescribed lift val- 
ue's. which arc very cle>se to the last computed case. 
For the 3 million point grid, the force coefficient val- 
ues at the prescribed lift conditions were interpolated 
directly from the 6 integer degree* eases within each 
drag polar. 
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Table 2 Definition of Required and Optional 
Cases for Drag Prediction Workshop 


Case 

Description 

Case 1 (Required) 

Mach = 0.73. C, = 0.500 

Single Point 


Case 2 (Required) 

Mach — 0.75 

Drag Polar 

n -3". 2”. I'M) 0 , r.2 1 ’ 

Cast' 3 (Optional) 

Mach = .50. .60. .70.. 75.. 76. .77.. 78. .80 

Constant Cj 

C, = 0.300 

Mat'll Sweep 


Case 1 (Optional) 

Mach = .30, .60, .70, .73. .76. .77. .78, .80 

Drag Rise Curves 

C{. - 0. 400. 0.300. 0.600, 


All cases were computed using the baseline grid 
(id) million points), and the medium grid (3 million 
points). Only the required cases were computed using 
the finest grid (13 million points) due to time con- 
straints. Table 3 depicts the results obtained for Case 
1 with the three different grids. The drag is soon to bo 
computed accurately by all Him' grids, although then' 
is a 10.6 count variation between the 3 grids. How- 
ever. the incidence at which the prescribed C\ ~ 0.5 is 
achieved is up to 0.6 degrees lower than that observed 
experimentally. This effect is more evident in the C\ 
versus incidence plot of Figure 6, where the computed 
lift values are consistently higher than the experimen- 
tal values. Since this discrepancy increases with the 
higher resolution grids, it cannot be attributed to a 
lack of grid resolution. The slope of the computed 
lift curve is about 5% higher than the experimentally 
determined slope, and is largely unaffected by grid res- 
olution. 

Figure 7 provides a comparison of computed surface 
pressure coefficients with experimental values at the 
experimentally proscribed Ci of 0.6 (where tin.' effects 
are more dramatic than at Cl — 0.5) as well as at the 
experimentally prescribed incidence of 0.93 degrees, at 
the 10.9 % span location. When the experimental inci- 
dence value is matched, the computed shock location 
is aft of the experimental values, and tin 4 computed 
lift is higher than the experimental value, while at the 
prescribed lift condition, the shock is further forward 
and the suction peak is lower than the experimental 
values. 

This bias in lift versus incidence was observed for a 
majority of the numerical solutions submitted to the 
workshop. 5 and thus might be attributed to a model 


geometry effect or a wind tunnel correction effect, al- 
though an exact cause has not been determined. When 
plotted as a drag polar. Cj versus Cp as shown in 
Figure 8, the results compare favorably with experi- 
mental data. Although the 4 drag polar was computed 
independently by both authors using the baseline grid, 
the results of both sets of computations were identical 
(as expected) and thus only one set of computations 
is shown for the baseline grid. The computational re- 
sults on this grid compare very well with experiment 
in the mid-range (near Cl = 0.5), while a slight over- 
prediction of drag is observed for low lift values, which 
decreases as the grid is refined. 

This behavior suggests an under-prediction of in- 
duced drag, possibly due to inadequate grid resolution 
in the tip region or elsewhere. The absolute drag lev- 
els have bet'll found to bo sensitive to the degree of 
grid refinement at the blunt trailing edge of the wing. 
The drag level is reduced by 1 counts when going from 
the 1.6 million point grid, which has 5 points on the 
trailing edge, to the 3 million point grid, which has 
17 points on the trailing edge. Internal studies by the 
second author using structured grids have shown that 
up to 33 points on the blunt trailing edge are required 
before the drag does not decrease any further. In the 
current grid generation environment, and without the 
aid of adaptive meshing techniques, the generation of 
highly refined trailing edge unstructured meshes has 
been found to be problematic, thus limiting our study 
in this area. 

Figure 9 provides an estimate of the induced drag 
factor, determined experimentally and computation- 
ally on the three meshes. 


Table 3 Results for Case 1; Experimental 
Values l:ONERA, 2:NLR, 3:DRA; GridV: Per- 
formed by first author, Gridl + : Performed by 
second author. Experimental data and 3 M 
point grid results are interpolated to specified 
Ci condition along drag polar. 


Case 

Ci, 

a 

Cp 

C\f 

Experiment 1 

0.3000 

+ .192* 

0.0289G 

-.1301 

Experiment 2 

0.5000 

f .153* 

0.02889 

-.1260 

Experiment 2 

0.5000 

+ .179* 

0.02793 

-.1371 

Grkll(l.6Mpts)' 

0.5004 

-.241” 

0.02921 

-.1549 

Gridl(l.6M pts U 

0.4995 

-,248 ( ' 

0.02899 

-.1548 

Grid2(3.0Mpta) 

0.5000 

-.417” 

0.02857 

-.1643 

Gri(i3(13Mpts) 

0.5003 

- .367” 

0.02813 

-.1657 
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Fig. 6 Comparison of Computed Lift as a 
function of Incidence for Three Different Grids 
versus Experimental Results 



Fig. 7 Comparison of Computed Surface Pres- 
sure Coefficients at Prescribed Lift and Pre- 
scribed Incidence versus Experimental Values 
for Baseline Grid at 40.9 % span location 



Fig. 8 Comparison of Computed versus Exper- 
imental Drag Polar for Mach=0.75 using Three 
Different Grids 



Fig. 9 Comparison of Computed versus Exper- 
imental Induced Drag Factor for TVIach~0.75 
using Three Different Grids 
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Fig. 10 Comparison of Computed ver- 
sus Experimental Idealized Profile Drag at 
Mach=0.75 using Three Different Grids 



Fig. 11 Comparison of Computed versus Ex- 
perimental Pitching Moment for Mach=0.75 
using Three Different Grids 

For Cr~ up to about 0.36. when the How is mostly 
attached, induced drag is underpredicted by approxi- 
mately 10 V. as determined by comparing the slopes 
of the computational and experimental curves (using 
a linear curve fit) in this region. Grid refinement ap- 
pears to have little' effect on the induced drag in this 
region. At tin' higher lift values, the 3 million point 
grid yields higher Cr. and lower Cp values, which is 
attributed to a slight delay in the amount of predicted 


How separation. Results for the 13 million point grid 
are not shown, since a fully converged solution could 
not be obtained at the highest incidence. On the other 
hand, it should be noted that the wind tunnel exper- 
iments used a boundary layer trip at 1 ~VX and 2oV 
chord on the upper and lower surfaces, while all cal- 
culations were performed in a fully turbulent mode. 
Examination of the generated eddy viscosity levels in 
the calculations reveals appreciable levels beginning 
between V/ to 7V chord. The exact influence of tran- 
sition location on overall computed force coefficients 
has not been quantified and requires further study. 

Figure 10 shows the idealized profile drag 19 which 
is defined by the formula: 

Cnr = Cn-C r ; 2 /(KAR) (1) 

where AR is the aspect ratio. Plotting Cpp generally 
results in a more compact representation of the data, 
allowing more expanded scales. It also highlights the 
characteristics at higher Cf, , where the drag polar be- 
comes non- parabolic duo to wave drag and separation. 
In the non-parabolic region, the error in drag is rela- 
tively large at a constant Cr. . 

The pitching moment is plotted as a function of Cr. 
in Figure 11 for all three grids versus experimental 
values. The pitching moment is substantially under- 
predicted with larger discrepancies observed for the re- 
fined grids. This is likely a result of the over-prediction 
of lift as a function of incidence, as mentioned earlier 
and illustrated in Figure 6. Because the computed 
shock location and suction peaks do not line up with 
experimental values, the predicted pitching moments 
can not be expected to be in good agreement with ex- 
perimental values. 

Figure 12 depicts the drag rise curve's obtained for 
Cases 3 and 4 on the baseline grid and the first refined 
grid (3 million points). Drag values are obtained at 
four different constant Cr, values for a range of Mach 
numbers. Drag value's art' predicted reasonably we'll 
except at the' highest lift and Mach number conditions. 
There appears to be no improvement in this area with 
increased grid resolution, which suggests issues such as 
transition and turbulence modeling may account for 
these discrepancies. However, since the two grids have 
comparable resolution in various areas of the domain, 
grid resolution issues still cannot be ruled out at this 
stage. 

The results obtained for Cases 3 and 1 can also be 
plotted at constant Mach number, as shown in the 
drag polar plots of Figure 13. The plots show simi- 
lar trends, with the drag being slightly overpredicted 
at low lift values on the coarser grid and with the re- 
fined grid achieving better agreement in these regions. 
For the higher Mach numbers, the drag is substantially 
underpredicted at the higher lift values. These discrep- 
ancies at the higher Mach numbers and lift conditions 
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point to an under-prediction of the extent of the sep- 
arated regions of flow in the numerical simulations. 
The comparison of idealized profile drag in Figure 10 
also suggests that the drag due to flow separation is 
not predicted accurately at the higher lift conditions. 
However, the character of the curves also suggest that 
the error may be due as well to the Cj offset (shown 
in Figure 6). Additional information concerning the 
regions of flow separation found in the wind tunnel 
would bo needed to more accurately quantify the na- 
ture of the error. 

The above results indicate that the current unstruc- 
tured mesh Navier-Stokes solver achieves a reasonably 
good predictive ability for tin' force coefficients on the 
baseline grid over the majority of the flow conditions 
considered. The overall agreement, particularly at the 
low lift values, is improved with added grid resolution, 
while the more extreme flow conditions which incur 
larger amounts of separation are more difficult to pre- 
dict accurately. On tin* other hand, tin* observed bias 
between computation and experiment in the lilt versus 
incidence values has an adverse affect on the prediction 
of pitching moment. While the source of this bias is 
not fully understood, it was observed for a majority of 
independent numerical simulations undertaken as part 
of the subject workshop 1 * and can likely be attributed 
to geometrical differences or wind tunnel corrections. 

The results presented in this paper involve a large 
number of individual steady-state cases. For example, 
on the baseline grid, a total of 72 individual cases were 
computed, as shown in Figure* 14. to enable the con- 
struction of Figures 8, 12, and 13. The majority of 
these cases were run from freest ream initial conditions 
for 500 multigrid cycles, while several cases particu- 
larly in the high Mach number and high lift regions 
were run 800 to 1000 cycle's to obtain fully converged 
results. The baseline cases (500 multigrid cycles) re- 
quired approximately 2.C hours of wall clock time on a 
cluster of 16 commodity PC processors. This enabled 
the entire set of 72 case's to be completed within a 
period of one week. This exercise illustrates the possi- 
bility of performing a large number of parameter runs, 
as is typically required in a design exercise, with a 
state-of-the-art unstructured solver on relatively inex- 
pensive parallel hardware'. 



Fig. 12 Comparison of Computed versus Ex- 
perimental Drag Rise Curves Three Different 
Cf. values on Two Different Grids 



Fig. 13 Comparison of Computed versus 
Experimental Drag Polars for Mach=0.6 and 
Mach— 0.8 on Two Different Grids 
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Fig. 14 Depiction of All 72 Individual Cases 
run on Baseline Grid Plotted in Drag Polar 
Format 


Conclusions 

A state-of-the-art unstructured multigrid Navier- 
Stokes solver lias demonstrated good drag predictive 
ability for a wing- body configuration in the transonic 
regime. Acceptable accuracy has been achieved on rel- 
atively coarse meshes of the order of several million 
grid points, while improved accuracy has been demon- 
strated with increased grid resolution. Grid resolution 
remains an important issue, and considerable exper- 
tise is required in specifying the distribution of grid 
resolution in order to achieve a good predictive ability 
without resorting to extremely large mesh sizes. These 
issues can be resolved to some degree by the use of au- 
tomatic grid adaptation procedures, which are planned 
for future work. The predictive ability of the numeri- 
cal scheme was found to degrade for How conditions 
involving larger amounts of How separation. Slight 
convergence degradation was observed on two of the 
grids for the cases involving increased How separation, 
white a fully converged result could not be obtained 
on the Hi lest grid (13 million points) for the highest 
lift case at a Mach number of 0.75. The current re- 
sults utilized the Spalart-Allmaras turbulence model 
exclusively, and the effect of other turbulence mod- 
els in this regime deserves additional consideration. 
The rapid convergence of the multigrid scheme cou- 
pled with the parallel implementation on commodity 
networked computer clusters has been shown to pro- 
duce a useful design tool with quick turnaround time. 
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